Quenches in quantum many-body systems: One-dimensional Bose-Hubbard model reexamined 
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When a quantum many-body system undergoes a quench, the time-averaged density-matrix p governs the 
time-averaged expectation value of any observable. It is therefore the key object to look at when comparing 
results with equilibrium predictions. We show that the weights of p can be efficiently computed with Lanc- 
zos diagonalization for relatively large Hilbert spaces. As an application, we investigate the crossover from 
perturbative to non-perturbative quenches in the nonintegrable Bose-Hubbard model: on finite systems, an ap- 
proximate Boltzmann distribution is observed for small quenches, while for larger ones the distributions do not 
follow standard equilibrium predictions. Studying thermodynamical features, such as the energy fluctuations 
and the entropy, show that p bears a memory of the initial state. 

PACS numbers: 05.70.Ln, 75.40.Mg, 67.85.Hj 
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Recent experiments d] in ultra-cold atoms have renewed 
the interest for the time-evolution of an isolated quantum 
many-body system after a sudden change of the Hamiltonian 
parameters, the so-called "quantum quench". Many questions 
arise from such a setup, among which are the relaxation to 
equilibrium statistics, the memory kept from the initial state, 
and the role of the integrability of the Hamiltonian. Analyt- 
ical and numerical results support different answers to these 
questions 10. HI B Htl, though most of them have shown that 
observables do not follow usual equilibrium predictions. As it 
has been pointed out |5, 6], looking at simple observables, yet 
experimentally accessible, might not be considered as suffi- 
cient to fully address these questions. Since time-evolution 
is unitary, there is no relaxation in the sense of a station- 
ary density-matrix, contrary to what can happen in a sub- 
system |6]. However, observables will fluctuate with time 
around some average. Standard definitions show that the time- 
averaged density-matrix p of the system governs any observ- 
able and its fluctuations. It is therefore deskable to have a 
systematic way of getting some information about p, and its 
associated thermodynamical-like quantities, in order to com- 
pare it with the density-matrices of equilibrium ensembles, 
such as the microcanonical or the canonical ensemble. 

In this paper, we show how Lanczos diagonalization (LD) 
enables one to calculate the weights of the time-averaged 
density-matrix. This method, which gives access to relatively 
large Hilbert spaces, is helpful when an analytical calcula- 
tion of the many-body wave-functions is lacking: this is, for 
instance, the case of nonintegrable models. As an applica- 
tion, the example of a quench in the one-dimensional Bose- 
Hubbard model (BHM) is revisited for the following reasons: 
(i) the model coiTesponds to realistic experiments [1], (ii) it 
is nonintegrable and it is usually believed that the redistribu- 
tion of momenta through scattering causes thermalization, (iii) 
complementary numerical results already exist |3], (iv) there 
is an equilibrium critical point demarcating two phases, and 
the latter can play a role in out-of-equilibrium physics. On 
finite systems, we show that there are two distinct regimes de- 
pending on the quench amplitude: in the perturbative regime, 
an approximate Boltzmann law is observed, while distribu- 
tions which do not belong to equilibrium ensembles emerge 



for large quenches. Moreover, we show that the mixed state 
p bears some memory of the initial state through its energy 
fluctuations and its entropy. 

We start by recalling |7] and introducing some definitions. 
From now on, the discussion will be restricted to finite-size 
systems of length L with no accidental degeneracy. We ad- 
dress the issue of the thermodynamical limit by looking at 
the scaling of observables with L, and by giving scaling ar- 
guments for the energy fluctuations. At time i < 0, the 
Hamiltonian is denoted by Tio and its eigenvectors and eigen- 
values by iV'n) and En- The system is prepared in some 
state I "00) J that usually is the ground-state of Hq. At i = 
0, the Hamiltonian is changed to 7i which eigenvalues and 
eigenvectors are cj„ and |0„). The time-evolving density- 
matrix of the whole system reads p{t) ~ X^nP" l*^") + 
E„<™ VP^[e"*"""*+'®"'" l-^") (0™|+/i.c.], with the rel- 
ative phases e„„i = 6'„ - e,-a, using 6'„ = Arg ((/)„|i/'o), and 
the frequencies ri„„j = lj^ - ^m- The p„ = |(-0o|?!'ri)P 
are the diagonal weights of the density-matrix, and they sat- 
isfy TlinPn ^ 1- As we are generally interested in the time- 
averaged expectation value of an observable O, we define 
O = limt^oo 7 /o Tr[p(s)0]ds = Y^nPnOnn, with the ma- 
trix elements Omn = (0n|O|(?!)m). Interestingly, averaging 
(0)o [with the notation {■)() — {t/jol ■ IV'o)] over random ini- 
tial phase differences 0„m gives back O, relating the time- 
averaging to the loss of information on the initial phases. Sim- 
ilarly, by averaging p{t) over time, one gets 
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which governs any time-averaged observable since O = 
Tr[pO]. Furthermore, it has been very recently shown |Q] 
that p is the experimentally relevant object to look at, and 
that the p„ weights enter in the microscopic expression of 
the work and heat done on the system in the quench. No- 
tice that the evolving state is a pure state so its von Neumann 
entropy S[p] = — Tr[plnp] is zero, while S\j5] is non-zero 
due to the loss of information induced by time-averaging. In 
addition, one must also look at the time-averaged fluctuations 
AO = [Tr[p(0-0)2]]V2 of the observables. We finally men- 
tion that, if O is diagonal in the |0„) basis, like the energy H, 
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the time-averaged expectations and fluctuations are fixed by 
the initial state: O = (O)o and AO = [((O - O f)oY^^. 

The difficulty for a given system is to compute the weights 
Pn or any expectation value. When there is no analytical ap- 
proach, as for the BHM, a possible solution is to resort to 
numerical techniques. In order to compute the p„, we notice 
that they enter in the expression of the (squared) fidelity |01 

Fit) = \A{t)\^ = 1 - 'iJ2n<mPnPmSm^[nn„^t/2]. This 

is the revival probability after a time t because we have 
A{t) = {ijj{t)\ipo), with \ip{t)) the time-evolving wave- 
function. A direct time-evolution calculation usually fails af- 
ter some time [3]. Our idea is to use spectral methods jgl fioll 
to get the Fourier transform Ajui) of the A{t) function. Con- 
trary to the approach of Ref. [l^ we notice that LD also 
gives a direct access to the Lehmann representation A{uj) = 
"^nPriH^ ~ + Eq) without a finite broadening, which 
induces an artificial decay of A{t). Hence, all the informa- 
tion we need to discuss the statistical features of p is included 
in A{uj), since both the energies and the weights are obtained. 
LD is not an exact method but is well adapted to low-energies, 
i.e. long times, and we give below a perturbative argument 
corroborating that the p„ have an overall decrease with ujn 
[see also ll 111 for cross-checking]. Hilbert spaces of sizes up 
to 10^ states will be studied in the following while our full 
diagonalizations are restricted to 5000 states. Lastly, spectral 
methods being much faster than time-evolution ones, one can 
scan a wide range of parameters. 

The short and long time behaviors of F{t) also contain in- 
formation about the p„ distribution |7]: at short times F{t) ~ 
1 — i^/r^ with r^^ = AE, the energy fluctuations. Phys- 
ically, the typical time r is the time after which the system 
has "escaped" from the initial state, and is the inverse of the 
centered width of A{uj). More generally, higher moments 
of die A{uj) function are defined by Mg = ([H - {n)o]'^)o, 
and are clearly fixed by the initial state. In practice, the mo- 
ments can also be independently computed with LD for q up 
to hundred by iteratively applying H on \tpo). The associ- 
ated sum rules are useful to cross-check the calculation of 
the spectrum. If one understands A{oj) as a probability dis- 
tribution, knowing all moments amounts to knowing the dis- 
tribution itself and would give back the exact p. This com- 
ment was put forward without proof in Ref. 4, together with 
a relevant discussion on the relation between these moments 
and generalized Gibbs ensembles. At long times, F{t) usu- 
ally fluctuates around its mean value = X^nPn 01- A 
qualitative interpretation of F is the "participation ratio" ^ 
that counts the number of eigenstates which contributes to 
time evolution. The typical fluctuations of the fidelity are 
(AF)^ 



F(i)2 - ^ 4 J2n<m PnPm- This quantity mea- 
sures the strength of the wavering of the evolving state be- 
tween getting back to ji/'o) or getting away from Ii/jq) ■ 

Qualitatively, a quench consists in projecting the initial 
state onto the spectrum of the Hamiltonian Ti. governing the 
dynamics. Straightforward results from perturbation theory in 
the quench amplitude illustrate the difference between small 
and large quenches: one expects a crossover between the two 
regimes. Writing 7i — Hq + XHi with A the quench am- 
plitude and Til the perturbing operator, the perturbed weights 



read, for A < 1, po - 1 - I]„^o ^"O, and pn^o - A^ft.„o, 
in which the notation hno = |(^/'n|'Wi|'0o)P/(£'n — -E-o)^ has 
been used. Meanwhile, the w„ are slightly shifted to order 
A and the eigenfunctions too. Thus, the p„ have an overall 
decrease with the excited energy and, increasing A induces 
a transfer of spectral weight from the "targeted" ground-state 
|(/)o) to other excited states. We get the scaling of several quan- 
tities to lowest order in A: Mq cx A^, 1 — F cx A^ and AF cx 
A^. As F > 0, these scalings will naturally fail for large A, 
signaling the crossover to the non-perturbative regime. In ad- 
dition, we mention that the mean-energy (E) is simply always 
linear in A, since we have (E) = {Ti)o = Eq + X{Hi)o- 

Application to a quench in the one- dimensional Bose- 
Hubbard model - We now study the BHM in a one- 
dimensional optical lattice which is a nonintegrable model: 
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fe,+fe]6,+i] + C//2^n,(n,-l). 



with 6| the operator creating a boson at site j and rij = b^jbj 
the local density. J is the kinetic energy scale while U is the 
magnitude of the onsite repulsion. In an optical lattice, the ra- 
tio U/ J can be tuned by changing the depth of the lattice and 
using Feshbach resonance yj. When the density of bosons is 
fixed at n = 1 and U is increased, the equilibrium phase dia- 
gram of the model displays a quantum phase transition from a 
superfluid phase to a Mott insulating phase in which particles 
are localized on each site. The critical point has been located 
at Uc — 3.3 J using numerics fl5l. The quenches are per- 
formed by changing the interaction parameter Ui — > Uf (we 
set J = 1 in the following), so we have A = {Uf — Ui)/2, 
and the perturbing operator Hi = J2j '^j ("-j ~ 1) is diagonal. 
Numerically, one must fix a maximum onsite occupancy. We 
take four as in Ref. 3 (for further details, see 1 1 1]). 

Since p features a mixed state, we call the {Ui, Uf) plane 
a state diagram. The Ui — Uf {X = 0) line splits this state 
diagram in two regions and the previous perturbative argu- 
ments should hold close to this line. The typical distributions 
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FIG. 1 : (Color online) Distributions of the p„ at four different points 
of the {Ui, Uf) state diagram. For the smallest size L — S, exact 
results are obtained by full diagonalization. 
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of the weights versus energy for four points of the state dia- 
gram are given in Fig. [1] two (a,c) with small quenches with 
parameters of the same (superfluid) equilibrium phase, and 
two (b,d) with large quenches, in which Uf "crosses" Uc in 
both directions. We observe that in the first two situations, 
for small A, the distributions are close to an exponential de- 
cay typical of a canonical ensemble. This result supports the 
evidence of a "thermalized" regime as found in Ref. 3, but 
on more general grounds since we directly have the distri- 
bution. Secondary peaks in Fig. [Tta) yield correction to this 
Boltzmann law. By looking at the cases of large quenches, we 
see that the distributions are strongly different from either the 
microcanonical or the canonical ensemble. When Uj — 20 
[Fig-Eb)], Mott excitations, corresponding to doubly occu- 
pied sites and roughly separated by U f, are clearly visible in 
the spectrum. Although the overall decay of the p„ is expo- 
nential, the distribution is very different from a Boltzmann 
law. This explains that many observables differ from the ones 
of an equilibrium system, and independently corroborates re- 
sults of Ref. 3. When Uf = 2 [Fig.ITJd)], the tai-geted spec- 
trum is nearly continuous and the distribution displays large 
weights around zero energy and a subexponential-like behav- 
ior [approximately exp(— (w„ — i?o)^) with 7 > 1]. This is 
again different from equilibrium predictions. The bump-like 
shape of the Uf — 2 distribution can be qualitatively under- 
stood from the fact that the ground-state energy increases with 
U in the BHM. As Eq > luq when Uf < Ui, the initial state is 
close in energy to some excited states of H and, according to 
the perturbative form of the p„, this favors their excitations by 
the quenching process. Another consequence is that the state 
diagram is expected to be non-symmetrical with respect to the 
Ui = Uf line. 

Crossover and finite size ejfects - To sketch the state dia- 
gram, maps of integrated quantities such as F, AF, and the 
entropy per particle s = S{p]/N are computed on a finite sys- 
tem with L ^ 12 and given in Fig. |2] The normalization of 
the entropy S [p] is motivated by the fact that we observe that 
it scales as N, plus some finite-size corrections. As suggested 
previously, observables display a crossover from the pertur- 
bative regime to a non-perturbative regime characterized by a 
significant enhancement of the weights of excited states. In 
order to evaluate the finite size effects on the crossover, we 
look at the scalings of F and AF for a cut along the Ui = 2 
line and increasing A. A first question is how the size of the 



perturbating regime evolves when increasing the length L. To 
address this question, we look at the evolution of two demar- 
cating points. One is associated with F and is inconclusive 
(for further details, see lHHl ). More interestingly, AF scales 
as A^ in the perturbative regime and the slope increases with 
L (see Fig. |3]l. At large A, AF is nearly flat and rapidly de- 
creases with L. In between, it passes through a maximum 
that defines a demarcating point Ac(F), and the correspond- 
ing AFc = AF(Ac(F)). The scaling of Xc{L) suggests a 
finite value in the thermodynamical limit [see Fig.|3jc)]. AFc 
can scale to a finite value but also to zero as a power-law [see 
Fig EIb)]. We notice that the latter situation would be in con- 
tradiction with a finite Ac and the fact that AF increases with 
L at low A. These results suggest that the perturbative regime 
survives in the thermodynamical hmit, but they remain ques- 
tionable. From the extrapolations, we find that the crossover 
survives for larger sizes, and could be experimentally relevant 
since experiments deal with finite systems. Notice that some 
of the numerics in Ref. 3 were done on larger systems. An- 
other question one can ask is the role of the critical point on 
the observed maximum of the fluctuations of the fidelity: one 
may define the "equilibrium expectation" Ac'' = {Uc — Ui)/2 
[resp. {Uf — Uc)/2] if one scans over Uf [resp. Ui] and com- 
pare it with the scalings of actual Ac. In Fig. [3tc) , the two 
are too close to be conclusive but for large Uij iflUl . the dif- 
ference is much substantial and Ac(F) even scales away from 
Ac"*. Thus, we infer that Uc certainly plays a role (see below 
and Fig. HI, but not on the crossover nor on the location of 
AFc. 

We now discuss some of the thermodynamical fea- 
tures of the mixed states described by the density-matrix 
p. Firstly, we ask whether the averaged energy is 
well-defined by looking at the relative energy fluctua- 
tions defined as AE/E = AF/((F) - Fq + AA^) = 

^JE^J{nln^j)o - {n'i)o{n^j)o/ E^i^Do to get rid of the ob- 
vious dependencies of (F) on Eq, N and A: what remains are 
the relative "squared density" fluctuations in the initial state. 
In the superfluid phase, we expect |13] the squared density- 
density correlations to have an asymptotic algebraic behavior 
(n|n|)o — (fT-f )o("-|)o ~ |« — il^", while they should be ex- 
ponential in the Mott phase e^'*^^'/^, with ^ the correlation 
length. On a chain of length L, we thus have AF = \\/Lg{L) 
with: (i) if a < 1, then g{L) - (ii) if a = 1, 
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FIG. 3: (Color online) (a) Cut along the Ui = 2 line showing a max- 
imum of AF between the perturbative and non-perturbative regimes 
of the quench, (b-c) Finite size scalings of AFc and Ac. See text for 
discussion. 




FIG. 4: (Color online) Rescaled relative energy fluctuations in the 
mixed state. They only depend on the features of the initial state. The 
slope of the curves vanishes close to the equilibriimi critical point Uc ■ 



g{L) ^ y/lnL and if a > 1 or ^ > 0, g{L) = const. As we 
have a > 1 in the superfluid phase of the ID BHM [13] and 
J2^{n'i)o - L, we find that AE/E = f {U, , L) / VL for any 
Ui. The f{Ui,L) function is computed with LD and plotted 
in Fig. m It shows a very good agreement with this scaling 



argument since / hardly depends on L. This scaling 
resembles the ones of the (micro)canonical ensembles but, 
one also notices that starting from a initial state with strong 
density fluctuations {a < 1) leads to anomalous scalings of 
the relative energy fluctuations. In the BHM, this could be 
achieved by introducing nearest neighbor repulsion lfl2ll . We 
also get the scaling of the typical time r ^ X~^L~^^^. This 
shows that, even if r scales to zero in the thermodynamical 
limit, it can be significantly long on large but finite systems 
for small A. More importantly, we find that two mixed states 
p can have the same energy (E) but with different A, since 
(E) = Eq + A(7ii)o. Hence, each of them originates from 
a different initial state and consequently, the two states have 
different energy fluctuations. Consequently, p keeps a mem- 
ory on the initial state. Another important thermodynami- 
cal feature is the entropy per particle s that continuously in- 
creases with A and reveals more significantly the underlying 
anisotropy of the state diagram [see Fig.|2l. We have checked 
that two mixed states with the same mean energy (E) have 
different entropies, so p also keeps a memory of the initial 
state through its entropy. From the very definition of the p„, 
this is not surprising. 

In conclusion, we have shown that the weights of the time- 
averaged density-matrix p can be obtained with LD. This pro- 
vides an observable-free description of the quench process, in 
particular for nonintegrable models. The method is applied to 
the ID BHM where it is shown that, on finite systems, there is 
a clear crossover from a perturbative regime, in which the dis- 
tribution is Boltzmann-like in the superfluid region, to distri- 
butions that are not predicted by equilibrium statistics ensem- 
bles. The state diagram has been mapped out in the {Ui, Uf) 
plane and finite size effects have been investigated. Lastly, we 
showed that the mixed state p has a well-defined energy and 
that it keeps a memory of the initial state through its energy 
fluctuations or its entropy. 

I thank T. Barthel, F. Heidrich-Meisner, T. Jolicoeur, D. 
Poilblanc and D. Ullmo for fruitful discussions. 
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ELECTRONIC PHYSICS AUXILIARY PUBLICATION SERVICE FOR: ON 
QUENCHES IN QUANTUM MANY-BODY SYSTEMS: THE 

ONE-DIMENSIONAL BOSE-HUBBARD MODEL REVISITED 

TYPICAL BEHAVIOR OF THE FIDELITY WITH TIME 
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FIG. 5: Typical behavior of the fidelity for a finite size system with L = 10 starting from Ui = 2 to Uf = 8 (A = 3). At short times: 
F{t) = 1 - t'^/r^ (here r = 0.14). For long times, F{t) fluctuates around its mean value F (here F = 0.135 and AF = 0.136). 



TECHNICAL DETAILS ON LANCZOS CALCULATIONS 



We use 200 Lanczos iterations to get the ground state and 1200 to get the Lehmann representation of A{uj). We do not use 
symmetries of the Hamiltonian except particle number conservation. With periodic boundary conditions, translational symme- 
tries induce some selection rules for the pn so their number is quite reduced. We have checked that Lanczos gives a good result 
by comparing it with exact results obtained by full diagonalization on a system with L = 8 (see Fig.|6] The largest Hilbert space 
size is 1331 1000 for L — lA for Lanczos diagonalization and 5475 for L = 8 for full diagonalization. Very similar results are 
obtained from systems with open boundary conditions. 




FIG. 6: Test on symmetries and effect of boundary conditions on the distribution of the pn. PBC stands for periodic boundary conditions while 
OBC is for open BC. 



ADDITIONAL RESULTS ON THE BOSE-HUBBARD MODEL 



Moments are related to the p„ and u;„ through Mq = J2nPn['-^n — {E)]''- They undergo a clear change of behavior with 
increasing A as shown in Fig.|7] 
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FIG. 7: First moments Mq (to the power of the A{ijS) function for a system with L = 10 and Ui = 2. There is a crossover from the 
perturbative result Mq ~ (gray lines) to a regime where Mq ~ A** at large A. 

We give below the behavior of F which goes from 1 when A = to small value when A is large. On a finite system, the 
second derivative (PF /d\^ crosses zero for a value \c{L) and we define the corresponding Ff. = F{Xc{L)). The scalings with 
1/L of these two quantities are given in Fig.[8] a linear scaling suggests that they are finite in the thermodynamical limit but 
power-law scalings going to zero also works for both, so studying this quantity is not very conclusive. Power-law scalings are 
however very slow, which means that even for large systems of length 100 or 1000 (experimentally relevant), the perturbative 
regime should survive. 

Notice that in the three other cuts of AF [Fig. [HI support that the same increase of AF with L in the perturbative regime as 
for the Ui ^ 2 case of the paper. 




FIG. 8: Left: Cut along the Ui — 2 axis of F. A linear scaling gives both a tinite value for Fc and Ac but a power-law one (going to zero) is 
also plausible. Right: Four cuts in the state diagram showing the scaling behavior of AF{L) as a function of A. The smallest size is i = 6 
and larger L — 13 except for Ui — 2 for which it is L = 14. Results for L = 6, 7, 8 are exact (full diagonalization) while larger sizes are 
obtained with Lanczos. The arrows indicate how AF increases or decreases with L. 



